%% SCRIPT TO PROCESS CFRADIAL DOW DATA FOR MARSHALL FIRE

close all
clear all


cd('~/Dropbox/MARSHALL_FIRE/')

%% Start by reading JSON file with fire perimtere data
fileName='marshall.json';
str = fileread(fileName); % dedicated for reading files as text
test=jsondecode(str);
figure(1);clf;
for aa=1:size(test.features.geometry.coordinates)
    figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot(pnow(:,1),pnow(:,2),'r')
end


%%


files=dir('./deployment1/cfrad*.nc');
allx=[];
ally=[];
allz=[];
alldbz=[];

xgrid=[(-5*10^4):50:2.5*10^4];
ygrid=[(-1.5*10^4):50:(.3*10^4)];
[xmt,ymt]=meshgrid(xgrid,ygrid);
[Vmaster,dbzmaster]=deal(nan(length(files),size(xmt,1),size(xmt,2)));
[evmaster]=deal(nan(length(files),1));
aaa=0;
evall=[]
for ff=1:length(files);

    fname=strcat(files(ff).folder,'/',files(ff).name);
    lats=ncread(fname,'latitude');
    lons=ncread(fname,'longitude');
    alt=ncread(fname,'altitude');
    vel=double(ncread(fname,'VEL'));
    velF=double(ncread(fname,'VEL_F'));
    dbzhc=double(ncread(fname,'DBZHC'));
    azimuths=double(ncread(fname,'azimuth'));
    elevation=double(ncread(fname,'elevation'));
    ranger=double(ncread(fname,'range'));
    width=double(ncread(fname,'WIDTH'));
    re=6.37*10^6;
    
    ts=double(ncread(fname,'time'));%time since volumne start
    starttime=char(ncread(fname,'time_coverage_start'));
    starttime=[starttime'];
    starttime(1,11)=' ';starttime(20:end)='';
    timer(ff)=datenum(starttime,'yyyy-mm-dd HH:MM:SS');

    clear xdist ydist dbznow

    for ev=1%:size(elevation)

        %%
        evnow=mean(elevation); % pull out the elevation for this scan.
        evall=cat(1,evall, evnow);
        if evnow>2 && evnow<4.3
            aaa=aaa+1;
            ae=(4/3).*(6.37*10^6); %earth radius
            z=(((ranger.^2) + (ae.^2) + (2*ranger.*ae.*(sind(evnow)))).^(1/2))-ae;
            hdist=  ae.*(asin((ranger.*cosd(evnow))./(ae+z)));%ae.*asin(rdist_HI.*(cosd(evnow))./(ae+z));
            [azmt,hdstmt]=meshgrid(azimuths,hdist); % create a mesh of azimuths and along ground distances
            [~,zdist]=meshgrid(azimuths,z);
            xdist=sind(azmt).*hdstmt; %x component of the horizontal distance
            ydist=cosd(azmt).*hdstmt;
            scalefactor=(re.*cosd(nanmean(lats))).*(pi./180);
            xdistlon=xdist./(scalefactor);
            ydistlat=ydist./(111180);

            %clean up the reflectivity data
            dbzhc=medfilt2(dbzhc);
            dbzhc(dbzhc<-5)=nan;
            dbzhc=medfilt2(dbzhc);
            width(isnan(dbzhc))=nan;

            %noise filter the velocity data
            dilnoise=imdilate(abs(del2(vel)),ones(5,5));
            mndbz=medfilt2(imerode(dbzhc,ones(5,5)));
            vel(isnan(dbzhc))=nan;

            figure(1);clf;
            subplot(1,2,1);cla;
            pcolor(xdist,ydist,width);shading flat; caxis([0 5])
            subplot(1,2,2);cla;
            pcolor(xdist,ydist,vel);shading flat; caxis([-10 20])
            
            %%
            addpath('~/Dropbox/MATLAB/SGP_LIDAR_CODES/')
            %             figure(1);clf;set(gcf,'pos',[100 100 1200 500],'color','w')
            %             sp1=subplot(1,2,1);cla;hold on;
            %             colormap(sp1,rbmapper(30,-10))
            %             pcolor(xdist,ydist,vel);shading flat;caxis([-10 30])
            %             cbh=colorbar;
            %             title(num2str(evnow),'fontsize',15)
            %             xlim([-1 3.5]*10^4);
            %             ylim([-2 1]*10^4);
            %             daspect([1 1 1])
            %             sp2=subplot(1,2,2);cla; hold on;
            %




            tic
            FV=TriScatteredInterp(xdist(:),ydist(:),vel(:));
            Vint=FV(xmt(:),ymt(:));

            FD=TriScatteredInterp(xdist(:),ydist(:),dbzhc(:));
            dbzint=FD(xmt(:),ymt(:));
            toc

            Vint=reshape(Vint,size(xmt));
            dbzint=reshape(dbzint,size(xmt));
            Vmaster(aaa,:,:)=Vint;
            dbzmaster(aaa,:,:)=dbzint;
            evmaster(aaa)=evnow;
            timemaster(aaa)=timer(ff);

            %             pcolor(xmt,ymt,(Vint));shading flat;caxis([-10 30])
            %             cbh=colorbar;
            %             colormap(sp2,rbmapper(30,-10))
            %             daspect([1 1 1])
            %             title(num2str(evnow),'fontsize',15)
            %             xlim([-1 3.5]*10^4);
            %             ylim([-2 1]*10^4);
            %             daspect([1 1 1])
            %             pause(.1)
            %             allx=[allx(:); xdist(~isnan(dbzhc))];
            %             ally=[ally(:); ydist(~isnan(dbzhc))];
            %             allz=[allz(:); zdist(~isnan(dbzhc))];
            %             alldbz=[alldbz(:); dbzhc(~isnan(dbzhc))];

        end
    end
end

%% Time mean velocity plot
evidx=find(evmaster>2 & evmaster<3.4);
% evidx=find(evmaster>3.4 & evmaster<5);
evidx=find(evmaster>2 & evmaster<5);

figure(22);clf;set(gcf,'color','w')
sp1=subplot(1,2,1);cla;hold on;
vbar=squeeze(nanmedian(Vmaster(evidx,:,:),1));
pcolor(xmt,ymt,vbar);shading flat;
caxis([-15 20])
cbh=colorbar;ylabel(cbh,'m s^{-1}','FontSize',15,'FontWeight','bold')
colormap(sp1,rbmapper(20,-15))
xlim([-1 4.5]*10^4);
ylim([-2.5 1]*10^4);
daspect([1 1 1])
hold on;
for aa=1:size(test.features.geometry.coordinates)
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot((pnow(:,1)-nanmean(lons)).*scalefactor,(pnow(:,2)-nanmean(lats)).*111180,'k','linewidth',2)
end
grid on;
box on;
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
xlabel('Dist [m]')
ylabel('Dist [m]')
xlim([-.5 1.5]*10^4);
ylim([-1 0]*10^4)
sp2=subplot(1,2,2);cla;hold on;
dbzx=squeeze(nanmax(dbzmaster(evidx,:,:),[],1));
pcolor(xmt,ymt,dbzx);shading flat;
caxis([0 30])
cbh=colorbar;
colormap(sp2,jet(12))
xlim([-1 3.5]*10^4);
ylim([-2 1]*10^4);
daspect([1 1 1])
hold on;
% figure(29)
for aa=1:size(test.features.geometry.coordinates)
    %     figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot((pnow(:,1)-nanmean(lons)).*scalefactor,(pnow(:,2)-nanmean(lats)).*111180,'k','linewidth',2)
end
grid on;
box on;
set(gca,'layer','top','linewidth',2)
xlabel('Dist [m]')
ylabel('Dist [m]')
xlim([-.5 1.5]*10^4);
ylim([-1 0]*10^4)
return
%% Examine temporal coefficient of variations for the wind data
Vcopy=Vmaster; Vcopy(abs(Vcopy)>30)=nan;
vstd=squeeze(nanstd(Vcopy(evidx,:,:),1));
vbar=squeeze(nanmean(Vcopy(evidx,:,:),1));
figure(80);clf;
sp1=subplot(1,2,1);cla; hold on;
pcolor(xmt,ymt,vbar);shading flat;
caxis([-15 20])
cbh=colorbar;ylabel(cbh,'m s^{-1}','FontSize',15,'FontWeight','bold')
colormap(sp1,rbmapper(20,-15))
xlim([-1 4.5]*10^4);
ylim([-2.5 1]*10^4);
daspect([1 1 1])
hold on;
% figure(29)
for aa=1:size(test.features.geometry.coordinates)
    %     figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot((pnow(:,1)-nanmean(lons)).*scalefactor,(pnow(:,2)-nanmean(lats)).*111180,'k','linewidth',2)
end
grid on;
box on;
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
xlabel('Dist [m]')
ylabel('Dist [m]')
xlim([-.5 1.85]*10^4);
ylim([-1.4 0]*10^4)

% sp2=subplot(2,2,2);cla;hold on;
% pcolor(xmt,ymt,vstd./abs(vbar));shading flat;
% caxis([-10 10])
% cbh=colorbar;ylabel(cbh,'m s^{-1}','FontSize',15,'FontWeight','bold')
% colormap(sp2,turbo)
% xlim([-1 4.5]*10^4);
% ylim([-2.5 1]*10^4);
% daspect([1 1 1])
%
% for aa=1:size(test.features.geometry.coordinates)
%     %     figure(1);hold on;
%     pnow=[squeeze(test.features.geometry.coordinates{aa})];
%     plot((pnow(:,1)-nanmean(lons)).*scalefactor,(pnow(:,2)-nanmean(lats)).*111180,'k','linewidth',2)
% end
% grid on;
% box on;
% set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
% xlabel('Dist [m]')
% ylabel('Dist [m]')
% xlim([-.5 1.85]*10^4);
% ylim([-1.4 0]*10^4)

%% Try the ratio of postive to negative occurences at a grid point
Vpositive=Vmaster(evidx,:,:);Vpositive(Vmaster(evidx,:,:)<0)=0; Vpositive(Vmaster(evidx,:,:)>0)=1;
Vnegative=Vmaster(evidx,:,:); Vnegative(Vmaster(evidx,:,:)>0)=0; Vnegative(Vmaster(evidx,:,:)<0)=1;
Vps=squeeze(nansum(Vpositive,1));
Vns=squeeze(nansum(Vnegative,1));

%% Ratios
% figure(20);clf;hold on;
sp2=subplot(1,2,2);cla; hold on;
pcolor(xmt,ymt,Vps./(Vps+Vns));shading flat;
for aa=1:size(test.features.geometry.coordinates)
    %     figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot((pnow(:,1)-nanmean(lons)).*scalefactor,(pnow(:,2)-nanmean(lats)).*111180,'k','linewidth',2)
end
grid on;
box on;
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
xlabel('Dist [m]')
ylabel('Dist [m]')
xlim([-.5 1.85]*10^4);
ylim([-1.4 0]*10^4)
title('fraction of time with positive radial wind')
cbh=colorbar;
colormap(sp2,parula(12))
daspect([1 1 1])


subplot(2,2,4);cla; hold on;
pcolor(xmt,ymt,Vns./(Vps+Vns));shading flat;
for aa=1:size(test.features.geometry.coordinates)
    %     figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot((pnow(:,1)-nanmean(lons)).*scalefactor,(pnow(:,2)-nanmean(lats)).*111180,'k','linewidth',2)
end
grid on;
box on;
title('fraction of time with negative radial wind')

set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
xlabel('Dist [m]')
ylabel('Dist [m]')
xlim([-.5 1.85]*10^4);
ylim([-1.4 0]*10^4)
return



%% Try making a google map version
addpath('~/Dropbox/MATLAB/')
Vcopy=Vmaster; Vcopy(abs(Vcopy)>30)=nan;
vstd=squeeze(nanstd(Vcopy(evidx,:,:),1));
vbar=squeeze(nanmean(Vcopy(evidx,:,:),1));
figure(80);clf;set(gcf,'color','w','pos',[100 100 1000 1000])
sp1=subplot(1,1,1);cla; hold on;
xlim(([-.5*10^4 1.85*10^4]./scalefactor)+lons(1));
ylim(([-1.4*10^4 0*10^4]./111180)+lats(1))
plot_google_map('MapType','hybrid');
pp=pcolor((xmt./scalefactor)+lons(1),(ymt./111180)+lats(1),vbar);shading flat;
% pp.FaceAlpha=.8;
% contour((xmt./scalefactor)+lons(1),(ymt./111180)+lats(1),vbar,[-30:5:30],'linewidth',2);shading flat;
caxis([-15 20])
cbh=colorbar;ylabel(cbh,'m s^{-1}','FontSize',15,'FontWeight','bold')
colormap(sp1,rbmapper(20,-15))

daspect([1 1 1])

hold on;
% figure(29)
for aa=1:size(test.features.geometry.coordinates)
    %     figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot((pnow(:,1)),(pnow(:,2)),'k','linewidth',2)
end
title(strcat(datestr(timemaster(evidx(1)),'HHMM'),'-',datestr((timemaster(evidx(end))),'HHMM'),{' UTC'}));
grid on;
box on;
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
xlabel('Dist [m]')
ylabel('Dist [m]')
xlim(([-.5*10^4 1.85*10^4]./scalefactor)+lons(1));
ylim(([-1.4*10^4 0*10^4]./111180)+lats(1))
%% Load station data and average over period of interest

mesos=dir('~/Dropbox/MARSHALL_FIRE/mesowest_data/*.csv');

for ms=1:length(mesos)

    T = readtable(strcat('/Users/nlareau/Dropbox/MARSHALL_FIRE/mesowest_data/',mesos(ms).name),'Delimiter',',');%,'Format','%s %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %s %f %f')
    jtime=datenum(T.Date_Time,'mm/dd/yyyy HH:MM');

    %just read header
    header=importdata(strcat('/Users/nlareau/Dropbox/MARSHALL_FIRE/mesowest_data/',mesos(ms).name),',',5)
    slat=header{3};slat=str2double(slat(12:end));
    slon=header{4};slon=str2double(slon(13:end));


    %compute the mean u,v wind components and plot as a vector on the map
    tidx=find(jtime>=timer(1) & jtime<=timer(end));
    uwind=-T.wind_speed_set_1.*sind(T.wind_direction_set_1);
    vwind=-T.wind_speed_set_1.*cosd(T.wind_direction_set_1);
    
    %%compute radial wind component
    azs=azimuth(lats(1),lons(1),slat,slon);
    vr=squeeze(uwind.*sind(azs))+ squeeze(vwind.*cosd(azs));


    %%
    ubar=nanmean(uwind(tidx));
    vbar=nanmean(vwind(tidx));
    vrbar=nanmean(vr(tidx))
    scalefactor=(re.*cosd(nanmean(lats))).*(pi./180);
    figure(80)
    %     qm=quiver(sp1,(slon-nanmean(lons)).*scalefactor,(slat-nanmean(lats)).*111180,ubar,vbar,200,'k')
    qm=quiver(sp1,(slon),(slat),ubar,vbar,.002,'k')

    sc=scatter(sp1,(slon),(slat),150,vrbar,'filled')
    sc.MarkerEdgeColor='k'
    sc.LineWidth=2;
    qm.LineWidth=3;
    qm.MaxHeadSize=3;
    
end

%%
%% Try making a google map version
addpath('~/Dropbox/MATLAB/')
Vcopy=Vmaster; Vcopy(abs(Vcopy)>30)=nan;
vstd=squeeze(nanstd(Vcopy(evidx,:,:),1));
vbar=squeeze(nanmean(Vcopy(evidx,:,:),1));
figure(81);clf;set(gcf,'color','w','pos',[100 100 1000 1000])
sp1=subplot(1,1,1);cla; hold on;
xlim(([-.5*10^4 1.85*10^4]./scalefactor)+lons(1));
ylim(([-1.4*10^4 0*10^4]./111180)+lats(1))
plot_google_map('MapType','hybrid');
pp=pcolor((xmt./scalefactor)+lons(1),(ymt./111180)+lats(1),Vps./(Vps+Vns));shading flat;
% pp.FaceAlpha=.8;
% contour((xmt./scalefactor)+lons(1),(ymt./111180)+lats(1),vbar,[-30:5:30],'linewidth',2);shading flat;
caxis([0 1])
cbh=colorbar;ylabel(cbh,'m s^{-1}','FontSize',15,'FontWeight','bold')
colormap(sp1,parula(12))

daspect([1 1 1])

hold on;
% figure(29)
for aa=1:size(test.features.geometry.coordinates)
    %     figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot((pnow(:,1)),(pnow(:,2)),'k','linewidth',2)
end
title(strcat(datestr(timemaster(evidx(1)),'HHMM'),'-',datestr((timemaster(evidx(end))),'HHMM'),{' UTC'}));
grid on;
box on;
set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
xlabel('Dist [m]')
ylabel('Dist [m]')
xlim(([-.5*10^4 1.85*10^4]./scalefactor)+lons(1));
ylim(([-1.4*10^4 0*10^4]./111180)+lats(1))

%% Load station and compute fraction of interval with positive vs. negative radial winds

mesos=dir('~/Dropbox/MARSHALL_FIRE/mesowest_data/*.csv');

for ms=1:length(mesos)

    T = readtable(strcat('/Users/nlareau/Dropbox/MARSHALL_FIRE/mesowest_data/',mesos(ms).name),'Delimiter',',');%,'Format','%s %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %s %f %f')
    jtime=datenum(T.Date_Time,'mm/dd/yyyy HH:MM');

    %just read header
    header=importdata(strcat('/Users/nlareau/Dropbox/MARSHALL_FIRE/mesowest_data/',mesos(ms).name),',',5)
    slat=header{3};slat=str2double(slat(12:end));
    slon=header{4};slon=str2double(slon(13:end));


    %compute the mean u,v wind components and plot as a vector on the map
    tidx=find(jtime>=timer(1) & jtime<=timer(end));
    uwind=-T.wind_speed_set_1.*sind(T.wind_direction_set_1);
    vwind=-T.wind_speed_set_1.*cosd(T.wind_direction_set_1);
    
    %%compute radial wind component
    azs=azimuth(lats(1),lons(1),slat,slon);
    vr=squeeze(uwind.*sind(azs))+ squeeze(vwind.*cosd(azs));
    vrcopy=vr(tidx);
    vrpos=vrcopy; vrpos(vrcopy<0)=0;vrpos(vrcopy>=0)=1;
    vrneg=vrcopy; vrneg(vrcopy>0)=0; vrneg(vrcopy<=0)=1;
    vrps=squeeze(nansum(vrpos));
    vrns=squeeze(nansum(vrneg));

    %%
    ubar=nanmean(uwind(tidx));
    vbar=nanmean(vwind(tidx));
    vrbar=nanmean(vr(tidx))
    scalefactor=(re.*cosd(nanmean(lats))).*(pi./180);
    figure(81)

    sc=scatter(sp1,(slon),(slat),150,vrps./(length(tidx)),'filled')
    sc.MarkerEdgeColor='k'
    sc.LineWidth=2;
    qm.LineWidth=3;
    qm.MaxHeadSize=3;
    
end